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Obtaining  accurate  computational  results  for  unsteady  low-Mach  number  flows  presents  significant 
challenges  for  numerical  schemes.  In  this  work,  we  present  a  novel  preconditioning  scheme  designed 
to  provide  good  conditioning  for  these  flows  for  all  Strouhal  numbers.  The  scheme  is  based  on  the 
convective-upwind  split-pressure  (CUSP)  scheme  which  naturally  facilitates  the  proper  scaling  of 
velocity  and  pressure  dissipation  terms  for  unsteady  low-Mach  number  flows.  Unlike  traditional 
matrix  dissipation  schemes,  the  preconditioned  CUSP  formulation  provides  scalar-like  efficiency  and 
simplicity.  Several  test  cases  are  presented  for  steady,  unsteady,  convection-dominated,  and  acoustic- 
dominated  flows  that  demonstrate  significant  improvements  in  accuracy  and  convergence  rate  over 
traditional  preconditioning  schemes  for  a  wide  range  of  flow  conditions. 


I.  Introduction 

Low  Mach  number  flows  present  significant  challenges  to  compressible  computational  fluid  dynamics  (CFD)  al¬ 
gorithms  in  terms  of  convergence  rate  and  accuracy.  Unlike  transonic  flows,  low  speed  flows  exhibit  widely  varying 
particle  and  acoustic  speeds,  which  have  a  negative  impact  on  convergence  rate.  Because  pseudo-time  scales  are  com¬ 
puted  based  on  the  fastest  (acoustic)  characteristic  speeds,  errors  propagating  at  much  lower  particle  velocities  take 
much  longer  to  convect  out  of  the  domain.  For  pseudo-time  approaches  that  rely  in  part  on  convection  of  errors  to 
reach  steady  state,  this  can  have  a  disastrous  effect  on  convergence  rate.  Even  if  steady  state  is  obtained  for  low  Mach 
number  flows  (after  many  iterations),  the  results  may  prove  inaccurate  due  to  improper  scaling  of  the  artificial  dissi¬ 
pation  operator.  This  has  been  observed  in  practice  and  in  theory  through  asymptotic  analysis  of  the  Navier- Stokes 
equations. 

Preconditioning  schemes  have  proven  highly  successful  at  improving  accuracy  and  convergence  for  steady  low 
Mach  number  flow  computations,  as  well  as  higher  Mach  number  flows  exhibiting  regions  of  low  speed  flow,  such 
as  boundary  layers  or  stagnation  points.  By  simultaneously  modifying  and  wave  speeds  in  pseudo-time  as  well  as 
artificial  dissipation  scaling,  preconditioning  can  extend  favorable  properties  observed  for  transonic  flows  to  low  Mach 
number  regimes.  Improvements  in  accuracy  and  convergence  can  be  dramatic,  as  demonstrated  by  many  researchers 
[1,  2,  3,  4].  Furthermore,  preconditioning  facilitates  the  use  of  a  variety  of  state  equations  necessary  for  real  gas  effects, 
multiphase  flows,  and  combustion  [5]. 

While  preconditioning  has  proven  highly  successful  for  steady  flows,  challenges  remain  to  extend  preconditioning 
to  unsteady  flows  for  which  propagation  of  vortical  and  acoustic  features  in  physical  time  becomes  a  complicating 
factor.  For  these  flows,  it  is  important  to  consider  the  effect  of  preconditioning  on  pressure  and  velocity  dissipation 
effects  separately.  At  low  Mach  numbers,  classical  “steady”  preconditioning  methods  increase  the  pressure  dissipation, 
while  reducing  the  velocity  and  temperature  field  dissipation,  thereby  insuring  well-conditioned  accuracy  at  all  speeds. 
The  situation  gets  murkier  for  unsteady  flows.  Here,  one  must  distinguish  between  vortex  propagation  and  acoustic 
wave  propagation.  For  the  former  case,  the  so-called  “low  Strouhal”  limit,  the  above  observation  remains  true  and  the 
steady  preconditioning  formulation  maintains  optimal  accuracy  at  all  speeds.  For  acoustic  propagation,  referred  to  as 
the  “high  Strouhal”  limit,  steady  preconditioning  results  in  overly  dissipative  pressure  dissipation,  while  the  velocity 
and  temperature  fields  are  optimally  handled.  On  the  other  hand,  when  no  preconditioning  is  used,  the  pressure  field 
is  optimally  handled,  while  the  velocity  and  temperature  fields  become  overly  dissipative. 
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In  order  to  address  the  dilemma  for  low-Mach  unsteady  flows  in  the  acoustic  limit,  blended  dissipation  schemes 
were  introduced  by  Xia  et  al.  [6]  for  the  scalar  dissipation  case  and  by  Potsdam  et  al.  [7]  for  the  matrix  dissipation 
case.  Further  work  on  unsteady  preconditioning  was  performed  by  Hosangadi  et  al.[8].  In  the  case  of  matrix  dissi¬ 
pation,  separate  treatment  of  pressure  and  velocity  becomes  difficult  due  to  the  simultaneous  coupling  of  these  terms. 
Nevertheless,  Potsdam  et  al.  [7]  demonstrated  that  is  it  still  possible. 

The  selection  of  different  scaling  for  the  pressure  and  velocity  diffusion  terms  is  perhaps  more  readily  accom¬ 
plished  via  the  AUSM  schemes  of  Liou  and  Steffen[9,  10,  1 1],  or  the  convective-upwind  split-pressure  (CUSP)  scheme 
of  Jameson[12,  13],  since  they  provide  distinct  formulations  for  these  terms.  For  the  steady  case,  Liou  provides  a 
modification  to  extend  AUSM  for  all  speeds,  which  he  terms  AUSM+-up[l  1].  Similarly,  Zha  et  al.[14,  15]  provides 
a  modification  to  CUSP  suitable  for  low  Mach  number  steady  flows,  although  it  appears  that  this  formulation  more 
closely  resembles  the  AUSM  framework. 

The  goal  of  this  work  is  to  explore  CUSP-based  preconditioning  schemes  suitable  for  all  speed  steady  and  un¬ 
steady  flows.  This  is  accomplished  via  the  addition  of  properly  scaled  pressure  diffusion  to  the  original  CUSP  scheme 
of  Jameson.  Flows  at  the  high  and  low  Strouhal  limits  are  explored  in  terms  of  accuracy  and  convergence.  Fur¬ 
thermore,  steady  cases  are  examined  to  verify  the  spatial  order  of  accuracy.  We  find  that  the  preconditioned  CUSP 
scheme  provides  accurate  results  and  good  convergence  at  all  flow  conditions  tested.  Moreover,  the  scheme  retains  the 
simplicity  and  economy  of  the  original  CUSP  formulation. 

The  paper  is  organized  as  follows:  First,  a  general  preconditioning  framework  is  outlined  based  on  primitive 
variables.  Modifications  for  steady  and  unsteady  flows  are  given.  Next,  the  formulation  of  the  new  preconditioned 
CUSP  scheme  is  given.  Following  this,  results  are  presented  for  steady  and  unsteady  flows  for  convective  and  acoustic 
dominated  cases  at  low  and  high  Strouhal  numbers.  Finally  preliminary  results  are  given  along  with  plans  for  future 
work. 


II.  General  Preconditioning  Framework 

The  purpose  of  this  section  is  to  formulate  a  general  framework  encompassing  most  preconditioning  schemes 
proposed  in  the  literature.  Preconditioning  schemes  have  been  proposed  by  many  researchers  to  improve  convergence 
and  accuracy  for  compressible  flows  using  a  pseudo-time  solution  procedure  [1,  2,  3,  4].  We  present  the  framework  in 
two-dimensions  for  simplicity,  but  extensions  to  three-dimensions  follow  readily.  In  subsequent  sections  we  use  this 
framework  to  describe  the  new  preconditioned  CUSP  scheme. 

Following  the  general  notation  of  Merkle  [5],  the  steady  inviscid  Euler  equations  are 
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Here,  Q  represents  the  vector  of  conserved  variables,  and  F  and  G  represent  the  inviscid  fluxes  in  the  x-  and  y- 
directions  respectively.  The  preconditioning  method  outlined  by  Merkle  involves  a  conversion  from  conserved  vari¬ 
ables  to  primitive  variables, 
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where  Qv  is  the  vector  of  primitive  variables.  The  use  of  primitive  variables  generally  simplifies  many  aspects  of  the 
present  derivation.  More  importantly,  primitive  variables  greatly  simplify  the  implementation  of  a  variety  of  equations 
of  state  and  real  gas  effects.  For  example,  enthalpy  is  often  tabulated  as  a  function  of  temperature  for  many  gases. 
With  conservative  variables  it  is  necessary  to  iteratively  find  temperature  from  enthalpy.  However,  since  temperature 
is  available  directly  from  the  primitive  variables,  the  enthalpy  can  be  extracted  with  little  effort  in  this  approach. 

Conversion  from  conserved  variables  to  primitive  variables  is  facilitated  via  the  Jacobian  matrix, 
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where  subscripts  p  and  T  denote  partial  differentiation  (e.g.  pp  =  dp/ dp).  Here  the  equation  of  state  is  formulated  as 

P  =  P(P,T),  h  =  h(p,T).  (5) 

Finite  volume  discretizations  of  Equation  1  involve  the  computation  of  directed  fluxes  at  volume  faces,  T  =  A  •  F, 
with  area  vector,  A  =  (Ax,  Ay)T .  The  wave  speeds  of  the  discrete  scheme  in  pseudo-time  are  the  eigenvalues  of  the 
product  of  T_1  and  the  directed  flux  Jacobian  with  respect  to  the  primitive  variables,  A  =  dT /dQv: 


X  (r  Aij  —  |  A  |  ( un ,  un ,  un  it  c) , 


(6) 


with  area  magnitude  \A\  =  y  A2  +  At*,  normal  velocity  un  =  u  •  n,  and  unit  normal,  n  =  A/|A|.  Here  c  is  the 
sound  speed,  defined  for  general  equations  of  state  as, 


c2  = 


phr 


phTpp  +  Pt(  1  -  php ) 


(7) 


Large  disparities  in  the  wave  speeds  of  Equation  6  result  in  poor  convergence.  The  goal  of  a  preconditioner  is  to 
reduce  the  wave  speed  disparities.  In  this  approach,  the  matrix  T,  arising  in  the  pseudo-time  derivative  is  replaced  by 
a  preconditioning  matrix,  Tp,  in  the  governing  equations: 


,9Q 


dF  dG_ 
dr  dx  +  dy 
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The  matrix  Tp  is  selected  such  that  the  eigenvalues  of  T  1 A  are  the  same  order  of  magnitude.  This  is  accomplished 
by  defining  Fp  as, 

p'p  0  0  pT 

up'p  P  0  upr 
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It  is  observed  that  Tp  is  identical  to  T,  except  that  a  new  term,  p'  replaces  pp.  In  terms  of  p'p,  the  eigenvalues  of  the 
preconditioned  system  become 

a  (r-U)  =  \a\ 

where 


U  n,  i  r/-r> ,  _ 


(9) 


(10) 


d=  phTpp  +  Pt(1  -  php) 
d!  =  phTp'p  + Pt{1  ~  php). 

The  selection  of  p'p  is  guided  by  the  principle  that  a  good  preconditioner  should  force  the  acoustic  wave  speeds  to 
be  the  same  order  of  magnitude  as  the  particle  speeds.  Following  Sankaran  [1]  and  Merkle  [5],  the  acoustic  waves  of 
the  preconditioned  system  in  Equation  involve  a  “preconditioned  sound  speed,”  Vp, 


V2  _  pfh 
p  ~  d! 


(ID 


If  one  selects  Vp  to  be  the  same  order  of  magnitude  as  the  particle  wave  speeds,  the  goal  of  preconditioning  is  achieved. 
Once  the  preconditioned  sound  speed  is  selected,  Equation  1 1  provides  a  definition  of  pp  in  terms  of  Vp : 


Pp 


V2 

Vp 


PA  1  ~  pK) 

phr 


(12) 


In  this  manner,  the  properties  of  the  preconditioner  depend  primarily  on  what  one  chooses  for  Vp.  The  following 
sections  provide  methods  to  determine  the  preconditioned  sound  speed  for  steady  and  unsteady  flows. 
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A.  Steady  Preconditioning 

For  steady  problems,  wave  information  propagates  at  convective  and  acoustic  speeds  in  pseudo-time  only.  In  this  case, 
the  preconditioned  sound  speed  is  set  to  the  magnitude  of  the  particle  wave  speed, 

Vp  =  vV+i’2.  (13) 

Here  the  superscript  “s”  denotes  that  this  definition  is  for  steady  flows.  At  low  speeds,  this  definition  of  V*  improves 
convergence  and  adjusts  the  scaling  of  the  artificial  dissipation  operator  [1].  For  supersonic  flows,  however,  precon¬ 
ditioning  is  unnecessary,  and  accuracy  actually  suffers  [1]  when  the  preconditioned  sound  speed  is  larger  than  the 
physical  sound  speed.  To  limit  the  preconditioned  sound  speed  for  high  speed  steady  flows,  we  choose, 

=  min  u 2  +  v 2,  (14) 

In  supersonic  regimes,  Equation  14  sets  the  preconditioned  sound  speed  to  the  physical  sound  speed.  This  reverts  the 
eigenvalues  to  their  non-preconditioned  form  and  turns  off  the  preconditioner. 


B.  Unsteady  Preconditioning 


For  unsteady  problems,  information  propagates  in  pseudo-time  and  physical  time.  As  the  relevant  physical  time  scale, 
At,  becomes  small,  physical  acoustic  waves  are  of  primary  consideration.  Hosangadi,  et  al.  [8]  showed  that  modi¬ 
fication  of  the  preconditioned  sound  speed  from  the  physical  value  in  this  case  can  cause  convergence  and  accuracy 
degradation.  In  this  discussion,  we  refer  to  flows  with  small  physical  time  scales,  At,  has  high  Strouhal  number  flows, 
where  the  Strouhal  number, 


Str  = 


L 

nAt\/u2  +  v2  ’ 


(15) 


represents  a  characteristic  non-dimensional  frequency.  High  Strouhal  number  flows  are  dominated  by  acoustic  wave 
action.  Flows  with  large  physical  time  scales  are  low  Strouhal  number  flows.  Taken  in  the  limit,  these  flows  revert  to 
steady  flows. 

In  this  light,  a  good  unsteady  preconditioner  should  turn  off  for  high  Strouhal  number  flows  to  allow  proper 
capturing  of  the  physical  acoustic  waves,  and  revert  to  the  steady  formulation  as  the  Strouhal  number  becomes  small. 
This  is  accomplished  by  defining  the  preconditioned  sound  speed  to  be 


V. pn  =  min 


max 


u2  +  u2,  Stry/ u 2  +  v2^j  ,  < 


(16) 


where  the  superscript  “u”  denotes  that  this  definition  is  for  steady  flows.  The  addition  of  the  Strouhal  number  in  this 
manner  allows  for  a  smooth  transition  between  steady  and  unsteady  (acoustically  driven)  flows. 

In  both  steady  and  unsteady  problems,  proper  treatment  of  stagnation  points  with  preconditioning  is  critical.  Stag¬ 
nation  points  are  common  in  external  aerodynamics  at  leading  and  trailing  edges  of  wings.  In  these  regions,  the  particle 
velocity  approaches  zero  and  the  definitions  of  Vp  and  Vp  become  ill-posed.  This  is  often  overcome  for  external  flows 
by  specifying  a  minimum  velocity  scale  proportional  to  the  free  stream  velocity.  In  stagnation  regions,  this  minimum 
velocity  scale  is  used  instead  of  the  particle  velocity  which  may  approach  zero.  For  other  types  of  flows,  such  as  inter¬ 
nal  flows,  a  free  stream  value  is  not  readily  available  and  using  the  particle  velocity  from  nearby  cells  as  a  reference 
often  works  well  in  practice.  For  the  cases  considered  in  this  work,  the  free  stream  specification  method  is  used. 


III.  Preconditioned  CUSP  Scheme 

The  previous  section  discussed  preconditioning  primarily  as  a  convergence  acceleration  tool.  The  other  critical 
aspect  of  preconditioning  is  its  effect  on  accuracy  through  the  scaling  of  numerical  dissipation  operators  needed  for 
stability  of  non-linear  systems  of  equations,  such  as  the  Euler  and  Navier-Stokes  equations.  This  can  be  seen  by 
examining  the  effect  of  preconditioning  on  a  finite  volume  numerical  flux  with  matrix  dissipation, 

^  (HQl)  +  F(Qr))  ~  r  p  iTp-^Ay  |  ( Qvr  ~  QvL )  •  (17) 

Extensive  research  has  been  performed  using  matrix  dissipation  with  both  steady  and  unsteady  preconditioners.  Hosan¬ 
gadi  et  al.  [8]  showed  through  asymptotic  analysis  the  general  behavior  for  these  schemes.  A  table  showing  the  effects 
of  Strouhal  number  and  Mach  number  on  the  accuracy  of  the  dissipation  for  the  pressure  and  velocity  fields  is  shown 
in  Table  1. 
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Steady  Low  Mach  Number  Limit 


Pressure  Field 

Velocity  Field 

No  Preconditioning 

O(M) 

0(1/M) 

Steady  Preconditioning 

0(1) 

0(1) 

Unsteady  Low  Mach/  Low-Str  Number  Limit 


Pressure  Field 

Velocity  Field 

No  Preconditioning 

O(M) 

0(1/M) 

Steady  Preconditioning 

0(1) 

0(1) 

Unsteady  Preconditioning 

0(1) 

0(1) 

Unsteady  Low  Mach/  High-Str  Number  Limit 


Pressure  Field 

Velocity  Field 

No  Preconditioning 

0(1) 

0(1/M) 

Steady  Preconditioning 

O(Str) 

0(1) 

Unsteady  Preconditioning 

0(1) 

0(1/M) 

Table  1.  Preconditioning  effects  on  artificial  dissipation  scaling  [8]. 


From  Table  1,  it  is  clear  that  matrix  dissipation  is  sufficient  for  steady  or  low  Strouhal  number  problems.  However, 
matrix  dissipation  accuracy  suffers  when  the  Strouhal  number  is  high  and  Mach  number  is  low.  In  these  cases,  steady 
preconditioning  is  capable  of  resolving  velocity  fields,  but  pressure  fields  suffer.  The  opposite  is  true  of  the  unsteady 
preconditioner.  To  remedy  the  situation,  Potsdam  [7]  created  a  “blended  matrix  scheme”  that  uses  selection  matrices 
to  combine  the  pressure  fields  of  the  unsteady  preconditioner  and  the  velocity  fields  of  the  steady  preconditioner. 

As  an  alternate  option,  we  consider  the  convective-upwind  split-pressure  (CUSP)  framework  proposed  by  Jame¬ 
son  [12].  The  CUSP  scheme  is  nearly  as  inexpensive  as  scalar  dissipation  but  retains  comparable  shock  capturing 
characteristics.  Furthermore,  it  naturally  splits  convective  and  pressure  effects,  which  proves  advantageous  in  a  pre¬ 
conditioning  context.  The  CUSP  scheme  is  given  by, 


D(Qr,  Ql )  =  ^a*|A|AQ  +  \f3  (T)Qr)  -  F(Ql)) 
where  A Q  =  Qr  —  Ql-  Equivalently  the  CUSP  scheme  may  be  expressed  as, 


(18) 


where 


D(Qr,  Ql)  =  ^a\A\cAQ  +  ^f3\A\QAun  +  ^/3ATP 


ac  =  a*c  +  /3un ,  Tp  — 


Axp 

AyP 

\A\unpt 


(19) 


(20) 


Here,  an  overbar  refers  to  the  arithmetic  average  of  right  and  left  states. 

In  the  CUSP  scheme,  a  and  (3  are  chosen  to  give  proper  up  winding  characteristics  in  subsonic  and  supersonic 
regimes,  while  maintaining  good  shock  capturing  for  transonic  regions: 


a  =  \Mn\ 

(21) 

r 

max(0, 2  Mn  —  1) 

if  0  <  Mn  <  1 

(3=1 

—  max(0,  2 Mn  +  1) 

if  -  1  <  Mn  <  0 

(22) 

\ 

sign(Mra) 

if  \Mn\  >  1 

Here,  Mn  =  un/c  is  the  local  directional  Mach  number. 

Unfortunately,  the  original  CUSP  scheme  above  encounters  similar  problems  as  the  matrix  scheme  at  low  Mach 
numbers.  In  order  to  determine  modifications  suitable  for  preconditioning,  we  note  that  Liou  et.  al.  [11]  observed 
similar  difficulties  with  AUSM.  To  solve  the  issue  in  AUSM,  Liou  introduced  two  new  forms  of  dissipation  into  the 
AUSM  scheme:  additional  pressure  dissipation  for  the  interfacial  Mach  number  and  additional  velocity  diffusion  to 
the  pressure  terms.  While  the  addition  of  pressure  dissipation  has  been  used  extensively  in  the  literature  [14,  15,  8] 
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with  much  success,  the  effects  of  velocity  dissipation  are  less  clear.  Since  the  focus  is  on  low  speed  flows,  velocity 
dissipation  is  omitted. 

Following  the  work  of  Liou,  additional  pressure  diffusion  may  be  added  to  the  CUSP  scheme  in  the  form 


Dp(Qr,  Ql)  =  7|  A\cQAp 


where 


7  = 


Kp  max(l  —  aMe  ,  0) 

fa 


pc 2 


fa  =  M0(2  -  M0),  Ml  =  min  (max  (Me  ,  Me  Str 2),  1), 

^2 


MP  = 


y/u2  +  v* 


Here,  Kv  is  taken  as  0.25  in  this  work.  The  final  modified  CUSP  scheme  is  then, 


D(Qr ,  Ql)  —  -<a|  A| cAQ  +  -f3\ A\QAun  +  -/3A T 0  +  Dp(Qr,  Ql) 


(23) 

(24) 

(25) 

(26) 

(27) 


While  CUSP  presents  advantages  in  terms  of  efficiency  over  the  matrix  scheme,  its  differentiation  does  not  natu¬ 
rally  lead  do  a  diagonally  dominant  left-hand- side  matrix,  which  causes  stability  problems  for  implicit  Gauss-Seidel 
schemes.  For  this  reason,  matrix  dissipation  is  always  used  on  the  left  hand  side  in  this  work. 


IV.  Results 

In  this  section,  we  present  steady  and  unsteady  results  to  assess  the  new  preconditioned  CUSP  scheme.  In  all 
cases  we  compare  with  established  matrix  dissipation  schemes.  First,  we  verify  that  all  preconditioning  schemes 
tested  have  been  implemented  properly  via  the  method  of  manufacture  solutions  (MMS).  A  steady  state  airfoil  is  also 
analyzed  for  validation  purposes.  Second,  in  order  to  assess  the  accuracy  for  low  Strouhal  and  low  Mach  number 
unsteady  flows,  we  compute  results  for  an  inviscid  propagating  vortex.  Finally,  flow  in  a  one-dimensional  pipe  with 
high-frequency  (in  time)  oscillating  back  pressure  is  computed.  The  high  pressure  oscillation  generates  acoustic  waves 
at  high  Strouhal  numbers.  The  results  from  Table  1  are  validated  using  matrix  dissipation  and  compared  to  the  new 
CUSP  algorithm  described  by  this  work.  Convergence  history  of  these  cases  is  provided  to  illustrate  the  effect  of 
low- speed  preconditioning  on  temporal  convergence. 

A.  Verification  and  Validation  for  Steady  Flows 

Verification  provides  a  robust  way  of  confirming  numeric  accuracy  and  code  implementation.  It  is  beneficial  to  use 
MMS  verification  any  time  a  new  scheme  is  being  introduced.  This  simple  check  was  performed  in  Fig.  1  using  steady 
state  MMS.  The  different  lines  and  symbols  found  throughout  Fig.  1  are  described  by  Table  2. 

The  MMS  verification  is  able  to  show  that  for  all  quantities  and  methods,  the  output  is  second  order  accurate 
and  behaves  properly.  It  is  interesting  to  observe  that  although  all  lines  are  second  order  accurate,  the  accuracy  is 
clearly  changing.  This  could  be  due  to  a  number  of  factors.  Regardless  of  whether  the  RHS  is  using  CUSP  or  matrix 
dissipation,  the  LHS  requires  the  use  of  the  diagonally  dominant  matrix  scheme.  For  the  parameters  of  source  terms 
chosen  to  define  the  MMS  solution,  the  optimum  preconditioned  pressure  additions  of  CUSP  could  be  providing  an 
excessive  amount  of  dissipation. 


Symbol 

Test  Scenario 

rn 

matrix  dissipation  with  No  preconditioning 

Rs 

matrix  dissipation  with  Steady  preconditioning 

CN 

CUSP  dissipation  with  No  preconditioning 

cs 

CUSP  dissipation  with  Steady  preconditioning 

Table  2.  List  of  test  scenarios  for  Fig  .1 
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RMS  Error  RMS  Error 


9 


Pressure 


10 


u  Velocity 


~20  40  60  80  100  120140 

1/h 


b)  u  Velocity 


d)  Temperature 


Figure  1.  Verification  plots  using  steady  state  MMS  for  primitive  variables  p ,  u,  v,  and  T  with  a  Mach  number  of  0.05.  Symbols  are  defined 
by  Table  2 
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To  get  general  validation  of  the  dissipation  schemes,  a  NACA  0012  airfoil  was  tested  at  a  relatively  low  Mach 
number.  Fig.  2  shows  the  pressure  contours  for  a  variety  of  cases.  A  free-stream  cut-off  velocity  was  specified 
for  preconditioning  at  the  stagnation  point  and  is  necessary  for  convergence.  Fig.  2c  shows  the  solution  without 
preconditioning.  The  stagnation  point  is  not  captured  in  the  detail  as  the  two  preconditioned  cases  found  in  Figs.  2a 
and  2b.  Near  the  stagnation  point,  the  particle  velocity  approaches  zero.  This  creates  a  very  large  disparity  in  the  wave 
speeds.  Although  the  solution  in  the  non-preconditioned  case  converges,  it  is  possible  for  the  overall  convergence 
criteria  to  be  reached  where  large  errors  are  still  propagating  locally  near  the  stagnation  point.  Local  preconditioning 
aids  in  overall  convergence,  as  well  as  the  convergence  at  the  stagnation  point. 

The  iterations  for  each  airfoil  is  also  reported  in  Fig.  2.  Without  preconditioning,  dramatic  efficiency  problems 
are  seen.  Implementation  of  CUSP  slows  down  convergence  slightly.  This  is  generally  attributed  to  the  difference 
between  the  LHS  and  RHS  in  the  implementation.  This  simple  case  highlights  how  preconditioning  is  able  to  increase 
the  accuracy  of  the  solution  as  well  as  aids  in  convergence. 


a)  matrix  with  Steady  Precon.  (345  Iterations)  b)  CUSP  with  Steady  Precon.  (584  Iterations) 


c)  matrix  without  Precon.  (2724  Iterations) 

Figure  2.  Pressure  contours  of  a  2D  NACA  0012  airfoil  with  free  stream  Mach  number  of  0.05  with  no  angle  of  attack. 


B.  Inviscid  Propagating  Vortex 

Inviscid  vortices  are  primarily  velocity  driven  problems.  Although  a  pressure  field  exists,  pressures  only  respond  to 
the  velocity  field.  An  inviscid  propagating  vortex  is  an  unsteady  problem,  where  some  vortex  travels  at  a  free  stream 
(E/oo)  velocity  in  time.  Ideally,  the  vortex  will  maintain  it’s  shape  with  time  and  simply  translate  through  the  domain. 
CFD  discretization  creates  diffusion  that  will  introduce  errors  into  the  solution.  To  test  this  phenomenon,  a  simple 
2D  vortex  is  used  to  initialize  the  domain.  Unsteady  time  steps  are  used  to  allow  the  vortex  to  propagate  across  the 
domain.  In  this  simplified  case,  it  is  assumed  the  free  stream  velocities  only  exist  in  the  x  direction. 

The  vortex  is  initialized  using  several  user  defined  parameters.  R  is  the  radius  of  the  vortex,  /?  is  the  strength  of 
the  vortex,  and  xc  and  yc  are  used  to  specify  it’s  starting  point.  Fig.  3  shows  a  sample  2D  initialization  of  the  vorticity 
over  a  square  domain.  The  following  equations  define  how  the  vortex  is  initialized  and  assume  an  ideal  gas  equation 
of  state. 
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Figure  3.  Initial  contour  of  vorticity  on  a  64  by  64  mesh. 


Su  =  -Uocp^e-7-2/2 

Sv  =  £/oo/3^e-r2/2  (28) 

ST  =  0.5(£/oo/3)2e-r7cp 

where 


uq  =  Uoo  +  u0  =  Sv,  T0  =  ^ 


P0  —  Poo 


Poo  — 


-Ra 


P()  —  P0  *  Pgas^ 0 


(29) 


y'(a;-xc)2  +  (y-7/c)2 

R 


Uoo  —  Moo  \/ T Pgas^oo 


Due  to  the  unsteady  nature  of  the  problem,  the  boundaries  of  the  domain  become  difficult  to  enforce.  The  exact 
solution  as  a  function  of  time  was  imposed  on  the  boundaries.  This  may  cause  issues  as  the  vortex  enters  or  leaves  the 
domain.  To  avoid  this  issue,  boundaries  were  placed  “far”  from  the  edges  of  the  vortex. 

In  the  absence  of  viscous  terms,  the  exact  solution  should  allow  the  vortex  to  travel  with  no  diffusion  at  the  free 
stream  velocity.  The  exact  solution  can  be  represented  by  the  initialization  equations  from  Eqs.  28  and  29  where  the 
center  is  moving.  The  new  center  is  found  as, 


^c(f)  —  ^c0  T  f/oot,  Pc(f)  —  UcO  (20) 

At  a  given  time,  the  exact  vorticity  can  be  found  to  be: 

UJ=E^e~r2/^2_r2)  (31) 

R 

For  unsteady  problems,  selection  of  the  physical  time  step  is  critical.  It  is  important  to  choose  a  time  step  to  capture 
the  flow  dynamics  efficiently  and  correctly.  The  “CFL”  number  relates  a  velocity,  time  step,  and  cell  width  to  help 
specify  the  time  step. 
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CFL  = 


Udt 

dx 


(32) 


where  U  is  some  velocity  scale,  dx  is  typically  taken  as  the  average  cell  width  and  dt  is  the  physical  time  step.  U  can 
be  chosen  to  be  either  the  particle  speed  ( CFLU )  or  the  acoustic  speed  ( CFLc ).  CFLU  =  1  allows  the  particles  to 
move  about  1  cell  per  time  step,  where  CFLc  =  1  allows  the  physical  acoustics  do  the  same.  These  two  situations  are 
of  current  interest  specifically  how  the  new  CUSP  modification  and  preconditioners  behave. 

In  order  to  validate  the  accuracy  of  Table  1,  as  well  as  the  implementation  of  the  preconditioners,  matrix  diffusion 
was  tested.  For  this  velocity  driven  case,  Table  1  indicates  that  using  no  preconditioner  and  the  unsteady  preconditioner 
in  the  high  Strouhal  limit  will  yield  reduced  accuracy. 

Fig.  4a  shows  the  data  for  a  “snapshot”  of  vorticity  at  some  time  after  the  vortex  is  allowed  to  translate  down 
stream  for  various  preconditioners.  The  vorticity  is  taken  at  the  centerline.  The  steady  preconditioner  captures  the 
vorticity  (or  velocity  field)  much  better  than  the  non-preconditioned  case.  The  introduction  of  the  Strouhal  number  in 
the  unsteady  case,  works  to  “blend”  the  steady  preconditioner  to  the  un-preconditioned  case.  Here,  the  time  step  is 
large  enough  that  the  high  Strouhal  limit  is  not  reached. 

Fig.  4b  shows  the  convergence  of  the  psuedo-time  terms  at  a  given  time  step.  The  unsteady  preconditioner  outper¬ 
forms  the  other  schemes.  By  converging  with  fewer  iterations,  the  unsteady  preconditioner  decreases  the  computation 
time  required.  Fig.  4b  highlights  just  one  of  many  time  steps.  As  the  number  of  time  steps  in  a  problem  grows,  this 
speed  up  becomes  increasingly  important.  The  results  from  Fig.  4  are  consistent  with  those  found  in  Hosangadi  et. 
al.  [8]  and  with  Table  1. 


Figure  4.  Plots  of  propagating  vortex  with  matrix  diffusion  on  a  square  domain  with  a  CFLU  =  1,  Str  =  20.4,  Mach  =  0.005,  and  20 
points  across  the  vortex. 

To  explore  the  effects  of  higher  Strouhal  numbers,  the  CFL  number  was  lowered.  Fig.  5  shows  a  similar  case  to 
Fig.  4,  but  with  a  CFLc  =  1  rather  than  CFLU  =  1.  The  reduction  in  the  time  step  causes  a  dramatic  increase  to  the 
Strouhal  number.  In  this  case,  the  unsteady  preconditioner  reduces  to  the  un-preconditioned  state,  consistent  with  its 
definition.  The  vorticity  plot  found  in  Fig.  5a  shows  this  behavior  and  highlights  the  unsteady  preconditioners  inability 
to  accurately  describe  the  velocity  field  at  high  Strouhal  numbers.  Again,  Fig.  5b  shows  that  the  steady  preconditioner 
suffers  to  efficiently  converge  in  this  highly  unsteady  case. 

The  unsteady  preconditioning  method  is  extremely  effective  at  converging  efficiently.  When  the  time  step  is 
very  small,  it  outperforms  the  steady  preconditioner.  As  time  steps  grow,  it  is  capable  of  reducing  to  the  steady 
preconditioner.  This  process  is  shown  by  Fig.  4b.  The  convergence  advantages  of  the  unsteady  preconditioner  are 
offset  in  this  scenario  by  the  drop  in  velocity  field  accuracy.  The  desire  to  capitalize  on  the  unsteady  preconditioners 
superior  convergence  while  not  loosing  accuracy  is  the  main  motivation  behind  the  implementation  of  a  preconditioned 
CUSP. 

In  order  to  showcase  the  effectiveness  of  CUSP,  the  cases  in  Figs.  4  and  5  were  re-run  using  the  newly  modified 
CUSP  diffusion  algorithm.  These  results  are  found  in  Figs.  6  and  7. 

Both  Figs.  6a  and  7a  show  accuracy  improvements  in  the  unsteady  preconditioned  case  when  compared  to  their 
matrix  counter  parts.  In  the  low  Mach  number  high  Strouhal  number  limit,  the  modified  CUSP  is  able  to  resolve 
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Figure  5.  Plots  of  propagating  vortex  with  matrix  diffusion  on  a  square  domain  with  a  CFLc  =  1,  Str  =  4076,  Mach  =  0.005,  and  20 
points  across  the  vortex. 


the  unsteady  issue  that  matrix  dissipation  exhibits.  Similar  to  the  airfoil  in  the  previous  section,  CUSP  shows  a 
slight  decrease  in  convergence  when  compared  with  matrix  dissipation.  This  is  most  likely  due  to  the  LHS  and  RHS 
discrepancies. 


Figure  6.  Plots  of  propagating  vortex  with  CUSP  diffusion  on  a  square  domain  with  a  CFLU  =  1,  Str  =  20.4,  Mach  =  0.005,  and  20 
points  across  the  vortex. 


C.  Oscillating  Back  Pressure  in  a  Pipe 

The  previous  case  highlighted  the  velocity  field.  In  order  to  isolate  the  pressure  field,  a  case  was  devised  to  purposefully 
introduce  pressure  waves  as  the  driving  factor.  A  simple  ID  pipe  was  modeled  using  a  2D  cell-centered  code  with 
uniform  cells  in  the  y-direction.  Pressure  on  the  outflow  was  oscillated  sinusoidally  in  time  to  create  pressure  waves. 
This  is  represented  by: 


Px=L  =  Poo(l  +  esin(wt))  (33) 

where  is  the  free-stream  pressure,  px=L  is  the  outlet  pressure,  w  is  the  frequency  that  the  pressure  is  driven,  t  is 
the  physical  time,  and  e  is  some  scale  factor. 

The  selection  of  e  was  set  so  the  overall  pressure  varied  less  than  1/2  of  the  dynamic  pressure.  This  prevents  the 
pressure  oscillation  from  reversing  the  flow  direction.  It  is  noteworthy  that  1/2  was  chosen  relatively  arbitrarily,  e  is 
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Figure  7.  Plots  of  propagating  vortex  with  CUSP  diffusion  on  a  square  domain  with  a  CFLc  =  1,  =  4076,  Mach  =  0.005,  and  20 

points  across  the  vortex. 


defined  below  by: 


e 


=  0.5 


0.5  Po 


(34) 


Another  dimensionless  parameter  is  introduced  that  relates  the  frequency  of  oscillation,  the  length  scale  of  the 
domain,  and  the  initial  flow  velocity.  The  length  scale  (length  of  the  pipe  L)  for  this  problem  was  set  to  one.  This 
parameter  (f 2 )  is  used  to  set  the  frequency  of  oscillation. 


^0 


(35) 


The  boundary  conditions  for  the  inlet  of  the  pipe  warrant  discussion.  It  is  necessary  to  fix  enough  of  the  state  to 
close  the  problem,  while  still  allowing  the  outflow  pressure  to  propagate.  One  method  is  to  use  non-reflecting  Riemann 
invariants.  Here,  stagnation  pressure  and  enthalpy  were  fixed  at  the  inlet.  These  allow  the  pressure  to  propagate  out  the 
inlet,  but  do  exhibit  some  reflective  behavior.  At  low  frequencies,  this  was  found  to  be  negligible.  At  high  frequencies, 
it  was  possible  to  observe  flow  behavior  before  the  information  arrived  at  the  inlet,  avoiding  the  issue  altogether. 

This  case  has  the  advantage  of  having  an  exact  solution  for  truly  incompressible  flows.  The  solutions  show  that  at 
a  given  time  step,  the  pressure  is  linear  in  space  while  the  velocity  is  constant.  The  exact  solution  is  presented  without 
derivation  below. 


Poouoo{l  +  U2) 


sin(wt )  —  Qcos(wt)  +  Ue  Uoot/L 


(36) 


x 

p(x,  t )  =  [esin{wt)  +  PooAxX^)]  ^  “  PooUoou(t)  (37) 

When  frequency  and  Mach  number  are  low,  the  incompressible  solution  should  be  a  good  approximation  of  the 
compressible  solution.  Figs.  8  and  9  show  this  situation  for  both  CUSP  and  matrix  dissipation.  Pressure  (at  the  center 
of  the  domain)  and  velocity  are  plotted  over  time.  Good  agreement  is  shown  between  the  exact  incompressible  solution 
for  all  methods,  again  validating  the  algorithms.  Figs.  8c  and  9c  reiterate  that  the  unsteady  preconditioner  is  clearly 
preferred  for  convergence  characteristics. 

To  understand  the  flow  behavior  as  frequency  increases,  it  is  necessary  to  explore  how  the  acoustic  wave  travel  in 
time  and  space.  The  period  of  the  outflow  pressure  oscillation  (Tw)  in  seconds  can  be  calculated  by: 


27 r  27 tL 

W  ^1A/[qqCqq 


(38) 


These  pressure  waves  travel  in  space  at  the  acoustic  speed  (c).  Multiplying  the  period  by  the  speed  at  which  the 
wave  travels,  yields  the  physical  period  ( Lw )  of  the  sinusoidal  pressure  wave. 
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Figure  8.  Plots  of  oscillating  back  pressure  with  matrix  diffusion  using  a  2D  mesh  of  128  by  4  with  a  CFLc  =  100,  Str  =  81.5, 
Mach  =  0.005,  ft  =  10. 
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Figure  9.  Plots  of  oscillating  back  pressure  with  CUSP  diffusion  using  a  2D  mesh  of  128  by  4  with  a  CFLc  =  100,  Str  =  81.5,  Mach 

0.005,  n  =  io. 
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Lw  -  Coo Tw  -  nM^  (39) 

When  frequencies  are  low  (i.e.  £2  =  10)  the  physical  period  of  the  wave  is  much  larger  than  the  domain  of 
the  problem.  In  these  cases,  the  pressure  will  appear  linear  and  the  incompressible  exact  solution  becomes  a  good 
approximation.  Lw  in  Figs.  8  and  9  was  approximately  125  times  larger  than  the  domain  length.  As  £2  increases 
the  solution  becomes  increasingly  invalid.  Fig.  10  shows  a  case  where  £2  =  100.  Here,  the  physical  period  is  about 
12.5  times  larger  than  the  domain  scale.  The  curvature  of  the  physically  acoustic  waves  is  now  no  longer  negligible. 
Fig.  10a  shows  the  pressures  at  the  center  of  the  domain.  It  is  clear  that  the  compressible  effect  of  the  pressure  waves 
is  beginning  to  manifest  in  the  solution. 


Figure  10.  Plots  of  oscillating  back  pressure  with  matrix  diffusion  using  a  2D  mesh  of  128  by  4  with  a  CFLc  =  1,  Str  =  8150,  Mach  = 

0.005,  n  =  too. 


To  explore  the  effects  of  the  acoustic  waves  over  the  domain  of  the  pipe,  £2  was  increased  to  4000.  This  creates 
a  physical  period  about  1/3  the  length  of  the  domain.  The  expected  solution  over  the  domain  is  a  perfectly  traced 
sinusoidal  pressure  wave  that  correlates  to  the  outflow  oscillations.  Figs.  11a  and  12a  show  the  pressure  waves,  as  they 
are  allowed  to  travel  leftwards  through  the  domain.  In  order  to  accurately  track  the  high  frequency  waves,  CFLc  was 
reduced  dramatically  to  0.025. 

Table  1  states  that  for  a  pressure  driven  problem,  the  steady  state  preconditioner  should  struggle  when  matrix 
dissipation  is  used.  This  is  consistent  with  the  wave  portrayed  in  Fig.  11a.  The  addition  of  the  modified  CUSP  scheme 
improves  the  steady  state  preconditioner  accuracy  in  this  regime. 

Even  though  CUSP  is  able  to  improve  the  accuracy  of  the  physical  pressure  waves  using  the  steady  preconditioner, 
the  steady  preconditioner  convergence  still  suffers.  This  is  highlighted  in  Figs,  lib  and  12b.  Unsteady  precondition¬ 
ing  in  combination  with  CUSP  is  capable  of  accurately  capturing  both  velocity  and  pressure  fields,  and  has  better 
convergence  characteristics. 


V.  Conclusions  and  Future  Work 

A  new  convective-upwind  split-pressure  (CUSP)  scheme  was  introduced  this  is  suitable  for  unsteady  flows  in 
the  low  Mach  number  high  Strouhal  number  regime.  The  new  scheme  follows  closely  the  original  CUSP  scheme 
of  Jameson,  but  with  additional  pressure  diffusion  to  improve  scaling  and  conditioning  at  low  Mach  numbers.  The 
behavior  of  the  preconditioned  CUSP  scheme  was  studied  for  acoustic  dominated  problems  by  setting  up  flow  in  a 
pipe  with  oscillating  back-pressure.  For  convection-dominated  flows,  the  scheme  was  evaluated  for  a  propagating 
inviscid  vortex.  The  new  CUSP  scheme  significantly  reduced  discretization  errors  in  both  cases.  It  was  also  shown 
that  the  introduction  of  the  Strouhal  number  into  the  preconditioning  definition  improves  efficiency  greatly  as  time 
steps  become  small.  Combining  the  new  CUSP  with  unsteady  preconditioning  is  robust  and  improves  accuracy  as 
well  as  efficiency  of  compressible  CFD  algorithms  for  low  speed  flows. 

Future  work  will  focus  on  the  linearization  of  the  CUSP  scheme  for  implicit  schemes  requiring  Jacobian  terms. 
All  results  is  this  work  were  obtained  using  matrix  dissipation  for  the  Jacobian  terms  regardless  of  the  right-hand  side 
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Figure  11.  Plots  of  oscillating  back  pressure  with  matrix  diffusion  using  a  2D  mesh  of  128  by  4  with  a  CFLc  =  0.025,  Str  =  326000, 
Mach  =  0.005,  Q  =  4000. 


Figure  12.  Plots  of  oscillating  back  pressure  with  CUSP  diffusion  using  a  2D  mesh  of  128  by  4  with  a  CFLc  =  0.025,  Str  =  326000, 
Mach  =  0.005,  Q  =  4000. 
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spatial  discretization.  This  strategy  was  adopted  to  preserve  diagonal  dominance.  However,  a  consistent  left-hand 
side  treatment  would  likely  further  improve  convergence.  Consistent  and  diagonally  dominant  linearization  of  the  new 
scheme  is  an  area  of  future  research.  In  addition,  future  work  will  focus  on  defining  the  Strouhal  number  based  on 
local  flow  conditions  rather  than  as  a  global  input,  which  requires  some  a  priori  knowledge  of  the  problem.  Local 
scaling  could  more  effectively  adapt  the  preconditioner  and  require  less  user  input. 
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